Collective Spin and Charge Excitations in Planar Aromatic Molecules 
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Employing high accuracy fixed node diffusion Monte Carlo (DMC) method we calculated the lowest triplet 
collective excitation (spin gap), as well as an upper bound for the singlet excitations (charge gap) in a series 
of charge neutral planar non-ladder aromatic compounds. Both excitation energies lie below the continuum of 
particle-hole excitation energies obtained from Hartree-Fock orbitals. Hence they can be interpreted as genuine 
bound states in the particle-hole channel. Assuming a resonating valence bond (RVB) ground state which has 
been recently suggested for sp^ bonded systems [ M. Marchi, et. al, Phys. Rev. Lett. 107, 086807 (2011)], 
offers a unified description of both excited states as two-spinon and doublon-holon bound states. We corroborate 
our interpretation, by Exact diagonalization study of a minimal model on finite honeycomb clusters. 



PACS numbers: 73.21.-b, 71.10.-w, 75.10.Kt 

Introduction: Correlation effects are characteristic of pi- 
TT conjugated systems composed essentially of hexagonal ar- 
rangements of sp'^ bonds 111]. Pauling took the initiative to 
describe the bonding in benzene (CeHe), a prototype of these 
systems, in terms of valence bonds (VB), focusing his atten- 
tion on the spin part of the bonding wave function, namely a 
singlet. Such singlet valence bonds can be formally described 
as the ground state of an effective Heisenberg exchange in- 
teraction JS1.S2, where J is the exchange integral between 
the overlapping atomic orbitals [2]. When the coordination 
number is low, in the above effective (model) Hamiltonian, 
the condition is more conducive for superposition of valence 
bond singlets to constitute the ground state - a unique oppor- 
tunity provided by three-fold coordination in these aromatic 
systems. Hence Pauling's formulation of the energy levels of 
molecules in terms of quantum mechanical superposition of 
valence bond configurations, the so called resonating valence 
bond (RVB) ||3|] becomes an important alternative rout to un- 
derstand energy levels of molecules. 



Some recent works ||4|-|6|] have combined the power of 
IVlonte Carlo methods with the basic notion of RVB [7], to 
construct variational wave functions in terms of geminals (de- 
terminants composed of "pairing" wave functions). Optimiza- 
tion of such RVB based many-body wave-functions, deter- 
mine the properties of sp'^ bonded systems with remarkable 
accuracy i^]^. ]V[ore specifically this technique captures the 
Kekule and Dewar contributions to the ground state of ben- 
zene 16] . Therefore, the notion of RVB in these systems is ca- 
pable of capturing interesting many-body effects in the ground 
state, at much lower computational cost, compared to more 
involved quantum chemical methods. The application of the 
above method to undoped graphene indicates that the ground 
state is a short range, gapped spin liquid i^ which agrees with 



other proposals based on the Hubbard model Isl-ilQl] . 

In view of the above mentioned evidences for a possible 
spin liquid ground state in planar sp'^ bonded systems, the 
next natural question would be about the nature of excited 
states. As a simple prototype which demonstrates the inad- 
equacy of single-particle description of energy states in this 
family of molecules, consider benzene, CqHq for which the 
IVlO would predict a singlet ground state for six pz electrons 
on the hexagonal ring. Within ]V[0 picture, the first singlet 
excited state (^i) and the first triplet excited state (Ti) are 
expected to be degenerate. However, observation of a remark- 
able splitting between the low-lying triplet and singlet excited 
states III2I] indicates the importance of correlation effects even 
in the excited states of these molecules. Since such correla- 
tion effects are based on local interactions, one expects the 
same picture to hold even in extremely extended members of 
this family, such as graphene U3\\ and carbon nano-tubes. The 
weak coupling (itinerant) limit of the graphene is well known 
to represent a Dirac liquid il6i], an d can be described by stan- 
dard si ngle -particle approach lll3ll . However, ab-initio calcu- 
lations [14j show that the strength of short-range part of the 
Coulomb interaction in these materials is ~ 10 eV, which is 
remarkably high and comparable to the estimated values of 
these parameters in conjugated polymers 11151] . For such large 
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values of Hubbard parameter U in these systems, emergence 
of a non Fermi liquid state, such as spin liquid [Sl |9|, [uJ] be- 
comes conceivable. 

In this paper, we investigate the nature of low-lying ex- 
cited states in small molecules belonging to the family of 
5p^ -bonded carbon systems. Here we employ the state of 
the art QIMC method to investigate the nature of many-body 
excitations in such hydrocarbons. This numerically accurate 
method suggests that the lowest excitation in such molecules 
is a triplet state, separated by a substantial gap from the next 
singlet excited state, for which we obtain an upper bound. We 
argue that these two lowest excited states, namely, Ti and Si 
can be naturally understood in terms of a picture based on 



spin-charge separation. This suggests that the ground state 
could be viewed as a resonating valence bond state, in agree- 
ment with a recent proposal by Marchi and coworkers |6]. 

Method: Considering computational cost and accuracy, 
Variational Monte Carlo (VMC) and Diffusion Monte Carlo 
(DMC) ||2Q|i algorithms are methods of choice for the cal- 
culation of many body properties of medium electronic sys- 
tems. These QMC methods can achieve chemical accuracy 
with a typical computational cost ranging from the second 
to fourth power of the number of particles ll24ll . In this pa- 
per we use these methods as implemented in CASINO pack- 
age to calculate spin and charge gap of some aromatic com- 
pounds. The CASINO code employs important sampling 
DMC method ||2ll l22|] to project out the many-body lowest 
energy state. In this method, the important sampled imagi- 
nary time Schrodinger equation is of the following form: 



/(R,^- 



Ar)= Jk{R, 



t^AT-R\t)f{R,t)dR\ (1) 



where /(R,t + Ar ) = ^^(R)V^(R,t + Ar ), ^t(R) is the 
trial wave function and V^(R, t + Ar ) is system wave func- 
tion. The kernel K{R, t + Ar ; R^ t) is the propagator. As 
Ar approaches to infinity, V^(R, t + Ar) tends to ground state 
in any sector corresponding to a definite set of quantum num- 
bers. For an efficient DMC calculation we need an optimized 
trial wave function. We used the multiplication of spin up and 
down Slater determinants and a Jastrow factor as a trial wave 
function: 



^, =^WDnri,...,r^)Z)J(ri,...,r^). 



(2) 



Here R = (ri,r2, ...rAr) denotes the spatial coordinates of 
all the electrons. The single-particle orbitals employed in 
the above Slater determinants have been constructed from 
Hartree-Fock (HF) mean field solutions which serve as a ref- 
erence basis for "free" particle-hole excitations. Note that this 
is not the exact Jastrow- Slater trial wave function form, as it is 
antisymmetric only with respect to the exchange of electrons 
with the same spin. Such wave functions can be used to obtain 
expectation values with lower computational cost for any spin 
independent operators ||2Q|i . CASINO uses Jastrow factors of 
the form proposed in Refs. 1221 123|]. We have taken into ac- 



count the electron-electron terms u, electron-nucleus terms x 
centered on the nuclei and 3 body electron-electron-nucleus 
terms/ in our calculations: 

A^-l A^ Anions A^ 

Anions A'- 1 N 

+ XI IZ IZ fi{ri,i,rjj,rij). (3) 
1=1 i=i j=i-\-i 

Optimization with respect to the parameters contained in the 
Jastrow factor was achieved by a VMC variance minimiza- 
tion procedure ||2QIi . After VMC optimization we used the so 
optimized wave function as a DMC trial wave function. Op- 
timization of Jastrow factors without optimizing orbitals did 



not affect the accuracy in our calculations. However, opti- 
mization of Jastrow factors provides a better trial wave func- 
tion for DMC calculation by making it more efficient. At the 
final stage of calculation, DMC projects out the ground state 
from this trial wave function. 

Using the above method, we calculate the many-body 
ground state in a given sector corresponding to the conserved 
total Sz and total number of particles TV. To extract informa- 
tion about spin-charge splitting from total energies, we pro- 
ceed as follows: Let Eo{N^^ N^) denote the ground state en- 
ergy for a system where No- is the number of electrons, each 
carrying spin ha/ 2 with a = ± corresponding to t and | spin 
orientations, respectively. TV = ^^ No- = N^ -\- N^, 2iS well 
as the total spin component, Sz = {N^ — N^)/2 are constants 
of motion and hence do not change the numerical projection 
by DMC procedure. Therefore quantum numbers {N^^^N^) 
appropriately label various sectors of the spectrum. 

Let us define the the spin gap (Ag) and charge gap (Ac) as. 

As = Eo{N^^l,N^-l)-Eo{N^,N^),, (4) 

-2Eo{N^,N^) (5) 

where (N^^N^) correspond to neutral system. In all com- 
pounds considered here, the total number of electrons, N, is 
even, so that the unpolarized configuration (i.e. the state with 
equal number of spin up and spin down electrons, N^ = N^) 
turns out to be the ground state. The energy Eq{N^^ N^) of 
this state can be calculated as follows: We generate a trial 
wave function from HF method with fixed total charge (neu- 
tral) and one spin multiplicity. Now to calculate the spin gap, 
Eq. dll) we flip one of the spins from, e.g. | sector, without 
altering the total charge. The lowest energy obtained by QMC 
procedure in this sector will correspond to ^o (^t + ^ ' ^i ~ ^ ) • 
Note that in this sector, the total charge is zero and spin mul- 
tiplicity is three. Note that since the spin and spatial sym- 
metries of the many-body Hamiltonian are not broken by 
HF solutions, the corresponding symmetry attributes are not 
changed by QMC projection. This means that the energy of 
the {N^ + 1, 7V| — 1) state will represent any of the three de- 
generate states belonging to the triplet representation of the 
SU{2) group. Spin gap defined above, represents the exact 
value of triplet excitation energy. 

Now let us discuss the physical meaning of the charge gap 
defined above: Imagine an infinitely large system, with equal 
number of t and ^ spin electrons. When an electron is moved 
from one point in the system to a distant location, the resulting 
excitation will be a doublon-holon pair. Ac is half of the aver- 
age energy of a pair, and hence can be interpreted as an upper 
bound for the energy of a single holon. To calculate the en- 
ergy of such doublon-holon pair, an approximate scheme is to 
isolate two small sub-systems surrounding the holon, and the 
doublon. In the absence of interactions, the doublon-holon 
energy will be given by the first term in the right hand side of 
Eq. ([5]). However, in reality there will be an attractive inter- 
action between them which lowers their true energy. There- 
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FIG. 1: Benzene CeHe, Phenanthrene C14H10, Pyrene CieHio, 
Coronene C24H12 



FIG. 2: Coronene with additional benzene ring C28H14 (benzo- 
coronene) 



TABLE I: Spin and charge excitations (eV) 



Compound Spin Gap 


Ac 


CeHe 


3.8(8) 


5.4(6) 


C14H10 


3.4(1) 


4.3(5) 


CieHio 


2.4(9) 


3.8(2) 


C24H12 


3.0(6) 


3.6(7) 


C28H14 


2.7(8) 


3.4(0) 



fore 2Ac defined above, is an upper bound for the energy of 
doublon-holon pair with respect to the neutral background. 
Because of the time-reversal symmetry of the Hamiltonian 
employed here, for such excitations based on charge fluctua- 
tions the spin orientation of the added/removed electron does 
not matter. 

Results: For five planar aromatic compounds depicted in 
Figs. lll2l we have calculated the above charge and spin gaps 
within the all electron fixed node DMC scheme. The re- 
sults are reported in Table [D The spin gaps obtained here are 
in good agreement with experimentally reported values lll2n . 
Also the charge gap Ac obtain here as an upper bound for the 
Si state agrees with existing results. For example in case of 
benzene, Ac = 5.4 eV represents a fair upper bound for the 
calculated result Es^ = 4.9 eV lVZ\. Therefore the method 
prescribed here to calculate the spin gap does indeed give the 
"lowest" excited state, and also the doublon-holon interpreta- 
tion employed here does represent a true upper bound for the 
energy of Si state. For geometry optimization as well as trial 
wave function generation, we used 6-31 IG** Gaussian basis 
which has been done by Gaussian 03 code 11911 . Note that 
we performed separate optimizations for ground, and excited 
states. All required energies are obtained with an accuracy 
better than ~ 5meV per atom. For each of the compounds 
reported in Table H and corresponding to each set of quantum 
numbers (A^^, A^^), we have optimized the geometry and the 
trial wave function constructed based on HF method. Then 
the Jastrow parameters have been optimized using variance 



minimization VMC. All reported DMC results are all-electron 
calculations, and we did not use any pseudo-potential. DMC 
time step is taken to be 0.002 Hartree"^. Optimized geome- 
tries have been verified to ensure they do not contain imagi- 
nary frequencies. 

To interpret the data in Table U let us represent them in a 
different way. In Fig. [3] we plot charge and spin gaps ver- 
sus the number of carbons. The new physical interpretation 
come about, when we also plot the tower of particle-hole ex- 
citations obtained from the (weakly correlated) Hartree-Fock 
theory for single particle states. This tower is the molecular 
analogue of the continuum of free particle-hole pairs. As can 
be seen by increasing the system size, the tower of particle- 
hole excitations approaches to a continuum. Moreover, the 
charge and spin gaps we obtain always remain below the the 
continuum of "free" particle-hole pairs. Therefore, they can 
be interpreted as the "bound state" of underlying free particle- 
hole pairs which are caused by many-body effects. First im- 
portant point which is suggested by this figure is that, a very 
large energy difference between lower edge of the tower, and 
the many-body states found here implies they are long-lived 
excitations which do not decay into the tower. Therefore, 
they can be associated with new quasi-particles. Note that 
the blue circles (Ac) is an upper bound for the true energy 
of a duoblon (holon), so that the true energy of the doublon 
state is even lower than Ac. The question is, what are these 
quasi-particles? 

Consider the lowest excited state Ti, which is a triplet 
many-body state for all system sizes considered here. The Ti 
state can be understood in terms of a simple RPA-like bound 
state formation in the triplet channel of particle-hole pairs. A 
short range repulsion of Hubbard type translates into the at- 
traction in the triplet particle-hole channel, and binds them 
together 1I10I1 . However the Si state whose exact location in 
Fig.[3]is somewhere between the blue circle and red square can 
not be understood in terms of simple RPA-like treatments, as 
the RPA in singlet channel predicts an anti-bound- state above 
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FIG. 3: Charge and spin gaps versus the number of carbon atoms in 
aromatic compounds studied here. To generate the "free" particle- 
hole continuum, we have used the Hartree-Fock orbitals. For all 
studied compounds, energy of spin excitation, Ti is below the sin- 
glet (charge) excitation 5*1 , and they both are below the continuum 
of free particle-hole excitations. 



the tower of free particle-hole states ilQII . Therefore the sec- 
ond excited state ^i is a genuine many-body effect, much be- 
yond the simple RPA like treatments. The method used here 
to obtain the upper bound for the singlet charge excitations 
suggests that the Si can be associated to an average energy of 
a doublon and a holon. To corroborate this claim further, let 
us use a simplified model Hamiltonian, which can capture the 
essence of the present QMC calculation in a more transparent 
way. First of all note that the minimal model which captures 
Ti state is a Hubbard model. Moreover, our earlier study of 
the particle-hole excitation spectrum in ID chains suggests 
that the singlet collective states below the particle-hole con- 
tinuum are controlled by the nearest neighbor Coulomb in- 
teraction 1 18]. Therefore the minimal effective model which 
captures both states is an extended Hubbard model. 



H = -t 



[hj)(^ 



4aCja+h.c.+/7 Y^ fij^rij^^V Y^ nifij. (6) 



Here z, j denote sites of a 2D honeycomb lattice and (i, j) 
indicates that they are nearest neighbors, ct creates an elec- 
tron in the Wannier state corresponding to the Pz orbital at 
site j. Here U and V denote the strength of on-site and near- 
est neighbor Coulomb interactions. Estimates of these param- 
eters based on the ab-initio methods indicates that even the 
screened of these parameters in graphene are substantial l\M 
and on the scale of [/ ~ 10 eV, which is comparable to corre- 
sponding estimates for smaller aromatic molecules f 15]. 

The result of the exact diagonalization for a 16-site hon- 
eycomb lattice is shown in Fig. IH The values of U = 9.3 
eV and t = 2.8 eV are adopted from Ref. 11141] . For the con- 
sidered range of V, the ground state always remains a total 
singlet state (So)- For small values of V, the first excited state 
is the Ti triplet, followed by a singlet excited state. Si . As 
V increases, the singlet excited state, ^i comes down and ap- 
proaches the energy of Ti excited state for F ^ 3 eV. Beyond 
this point, the first excited state will be a singlet state. Think- 
ing from the limit of very large molecules, the Si state will 




FIG. 4: (Color online) The excitation energy for a 16-site cluster 
corresponding to CieHio in the extended Hubbard model. Values of 
Hubbard U and t taken from Ref. iil4i1 are in eV. The ground state 
always remains a singlet. The first excited states for small values of 
V are triplet. By increasing V, the order of Ti (red) and Si (blue) 
is switched, and beyond V ^ 3 eV, Si will be the first excited state. 
All energies in this figure are in eV. 



have no analogue in terms of plasmon oscillations. Since, first 
of all, plasmon oscillations require non-neutral system OSll 
Secondly, long range Coulomb interaction makes the singlet 
branch either an acoustic plasmon (for 2D coulomb repulsion) 
or a gapped pi-plasmon (3D coulomb repulsion) branch. So 
the Si state can not be interpreted as molecular analogue of 
plasmon mode. On the other hand, the decreasing behavior of 
the Si energy with V is consistent with a doublon-holon inter- 
pretation: The repulsion V among the electrons will become 
attraction —V between the doublon-holon pair, and increasing 
V will lower their energy. 

Summary and discussions: We have used ab-initio QMC 
method to obtain an accurate excited state Ti and an upper 
bound for the Ti . We then used exact diagonalization to study 
a minimal model which captures the same set of excitations. 
Assuming RVB ground state |t6|], offers a unified understand- 
ing of both states. In this scenario, the Ti can be understood as 
the energy required to break a singlet in the RVB background 
and render it triplet lIlTII . Moreover, the Si can be attributed to 
a holon-doublon pair created by removing one electron from 
one carbon site, and placing it in the Pz orbital of another car- 
bon site. Such charge fluctuations are allowed because the 
on-site Coulomb energy U is finite. In this picture, the de- 
crease in ^i energy by increase in V becomes quite natural. 
This interpretation can be a possible description of the col- 
lective char ge e xcitations observed in thick multi-wall carbon 
nano-tubes 11260 . 
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